/*
 * KineticProblem.cpp
 *
 *  Created on: 18 Oct 2011
 *      Author: allan
 */

#include "KineticProblem.h"

// GeoReact includes
#include <Core/KineticSolver.h>

KineticProblem::KineticProblem(KineticSolver& kinetic_solver, double T, double P, const VectorXd& n) :
kinetic_solver(kinetic_solver), T(T), P(P), n(n)
{}

void KineticProblem::Function(double t, const VectorXd& ujnk, VectorXd& Rj_Rk)
{
	// Auxiliary variables
	const unsigned Nj = kinetic_solver.cannonical_system.GetNumPrimarySpecies();
	const unsigned Nk = kinetic_solver.cannonical_system.GetNumKineticSpecies();
	const Multiphase& multiphase = kinetic_solver.cannonical_system.GetMultiphase();
	const ReactionSystem& kinetic_reactions = kinetic_solver.cannonical_system.GetKineticReactions();
	
	// Solve the total equilibrium problem
	kinetic_solver.TotalEquilibrate(ujnk, T, P, n);
	
	// Calculate the activities of the species (a)
	const VectorXd a = multiphase.Activities(T, P, n);
	
	// Calculate the kinetic rates of kinetic reactions (rk)
	const VectorXd rk = kinetic_reactions.KineticRates(T, P, n, a);
	
	// Calculate the kinetic rates of the primary and kinetic species (Rj_Rk)
	Rj_Rk.segment(00, Nj) = kinetic_solver.cannonical_system.KineticRatesPrimarySpecies(rk);
	Rj_Rk.segment(Nj, Nk) = kinetic_solver.cannonical_system.KineticRatesKineticSpecies(rk);
}
